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SUMMARY 


Three-dimensional unsteady viscous effects can significantly Influence the 
performance of fixed and rotary wing aircraft. These effects are important in 
both flows about helicopter rotors in forward flight and flows about three- 
dimensional (swept and tapered) supercritical wings. A computational procedure 
for calculating such flow field is developed, and therefore would be of great 
value in the design process as well as in understanding the corresponding flow 

phenomena . 
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INTRODUCTION 


In recent years there has been increased attention given to three-dimensional 
unsteady aerodynamics. Such flows manifest themselves over fixed wing and rotary 
wing aircraft. In regard to rotary wing aircraft, the helicopter operates in an 
unsteady environment. The flow about a helicopter rotor in forward flight is 
periodic as the blade passes through the rotor disc and the flow is characterized 
by its unsteady three-dimensional nature. Further, unsteady effects result from 
the wake vortex interaction due to the shed vortex of the preceding blade passing 
in the vicinity of the subject blade. In regard to fixed wing aircraft they are 
designed to be nominally steady, unsteady effects are introduced through either 
control surface motions or induced external oscillations. For both wing types 
these phenomena are observed throughout the Mach number regime; subsonic, 
transonic, and supersonic. The unsteady effects will influence loss levels as 
well as lift and moment coefficients which in turn influence the aeroelastic wing 
loading. For both type of wings the near wing flow may contain transonic shock 
wave boundary layer interactions as well as significant regions of reversed flow 
in the streamwise and spanwise directions. Obviously, a viscous analysis capable 
of treating unsteady, three-dimensional flows which may contain shock wave 
boundary layer interactions as well as regions of spanwise and streamwise reversed 
flow would be a significant aid to both the design and research engineer. 

Concurrent with these needs there has been an increased effort to better 
understand these phenomena by conducting three-dimensional unsteady wing 
experimental programs (cf. Refs. 1, 2) and applying inviscid computational 
procedures to these programs (e.g. Refs. 3, 4). The results of the numerical 
calculations in conjunction with the experimental data indicate that the observed 
phenomena are strongly influenced by viscous effects near the body surface which 
are not accounted for by the inviscid predictions. These viscous effects are 
concentrated within a region that is predominantly thin except for localized 
regions of reverse flow in the streamwise and/or spanwise directions. Hence, 
there is clearly a need to compute these viscous effects in an efficient and 

economical manner . 

There are several possible approaches available for computing 
three-dimensional viscous flows, ranging from empirical models to sophisticated 
treatments based on the solution of the three-dimensional time -dependent 
Navier-Stokes equations. Due to the complex structure of the flow the empirical 
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approach is too restrictive. At the other end of the spectrum is the 
three-dimensional time -dependent Navier- Stokes analysis. Even though such 
procedures have been developed at SRA and have been applied successfully to a 
variety of problems (e.g. Ref. 5), such a technique is not required for many of 
the viscous layer type problems occurring on wings in which the static pressure is 
sensibly constant across the viscous layer. Therefore, an approach is sought 
which allows the static pressure to be imposed at the boundary layer edge, but 
which can be used in three-dimensional flows having streamwise and/or crossflow 
separation. 

A new computational procedure specifically designed to compute flow fields in 
which streamwise and/or spanwise separation is present, but in which the pressure 
is sensibly constant across the boundary layer has recently been developed 
(Ref. 6). This bridges the gap between the inviscid/boundary layer and 
Navier-Stokes approaches in that it is of sufficient generality to compute regions 
of reverse flow yet due to the imposition of pressure is considerably more 
economical than a full three-dimensional Navier-Stokes procedure. This technique 
was adopted to treat three-dimensional unsteady turbulent flows (Ref. 6) . The 
results of this study are now described briefly. 

Ref. 6 describes and demonstrates the implementation of a computer code for 
the efficient solution of three-dimensional time -dependent viscous flows on fixed 
and rotary wing aircraft. The numerical technique used is the Linearized Block 
Implicit (LBI) technique of Briley and McDonald (Refs. 7) in conjunction with QR 
operator technique (Refs. 8 and 9). This combination numerically solves the 
present approximate form of the turbulent Navier-Stokes equations which are 
derived for nonorthogonal coordinates in generalized tensor form. The rationale 
for the choice of this approach is discussed in detail in Refs. 6, 8 and 9. 

The basic assumption made in the derivation of the governing equations is that 
the pressure does not vary normal to the shear layer, and is obtained from an 
inviscid analysis. Inherent in this assumption is that the shear layer is thin. 
Generally speaking, the boundary layer remains thin unless catastrophic flow 
separation occurs or the flow at the wing or rotor tip is considered. However, 
the analysis would apply to most of the wing or rotor under a range of operating 
conditions and thus represents an important tool. 

It is also assumed at present that the stagnation temperature, T 0 , is 
constant. This assumption is a good approximation for the flow fields considered 
as discussed in Ref. 5, and is included here only for purposes of computer run 
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economy. The full energy equation could equally well have been used in the 
analysis with consequent increase in computer run time. When the total 
temperature is assumed constant, the equation of state relates the density p to 
the velocity components u and w by an algebraic equation. The resulting 
formulation involves only the three velocity components, u, w and v and three 
equations, the streamwise and spanwise momentum equations and the continuity 
equation. Hence, a block- three system is considered. If T 0 were calculated via 
the energy equation, a block- four system would result due to the inclusion of the 
temperature as an additional unknown and thus would result in an increase in 
computer run time. 

For turbulent flows, a two -layer mixing length model is employed and its 
formulation in generalized tensor notation is given. A novel method is employed 
for solving the continuity equation in conjunction with the momentum equations. 

In Ref. 6, a complete description of the computational procedure is given, 
including coordinate systems, governing equations, turbulence model, and numerical 
technique, i.e. QR operator and Linearized Block Implicit schemes. The general 
outline of the computer code is also described. 

The computational procedure has been validated by conducting computations with 
the numerical method referred to above and comparing it to the experimental data 
of Karlsson (Ref. 10), i.e. two-dimensional unsteady oscillating turbulent flow 
over a flat plate. Two-dimensional calculations were performed and the results 
agree both qualitatively and quantitatively with the data. Thereafter, the 
analogous three-dimensional case was considered which was obtained by a coordinate 
rotation to yield the flow over a plate skewed at 45° to the frees tream direction. 
The results of this computation also agree well with both the two-dimensional 
results and Karlsson' s data, hence validating the computational procedure in three 
dimensions. In addition, new inflow boundary conditions were developed and an 
explanation was proposed to resolve the controversy concerning other previously 
reported predictions of the skin friction phase lead angle as a function of 
reduced frequency. These results are described in detail in Ref. 6. 

In this report, a description is given of the extension of the methodology to 
treat two-dimensional airfoils and three-dimensional wings. In order to achieve 
these applications a generalized three-dimensional nonorthogonal geometry 
formulation was developed. Also, the procedure was modified to calculate the flow 
through the wake. Examples of calculations conducted are presented. 
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DISCUSSION 


To date, a computational procedure for three-dimensional unsteady viscous 
flows has been validated for three-dimensional unsteady oscillatory flow over a 
planar surface. The method was shown to be efficient and gave both qualitative 
and quantitative agreement with the experimental data. The major goal of the 
current effort was to extend this procedure so that it may be used routinely as an 
aid in the design of realistic fixed and rotary wing aircraft. This entailed in 
part the validation and extension of options in the computer code that were not 
exercised as yet and also required the incorporation of additional capabilities 
that would allow one to consider a more general class of flow phenomena consistent 
with current interests to both NASA and industry. 

A major task required to complete the long term goal was to extend the 
geometrical capability of the computer code. The items considered were validation 
of the existing generalized nonorthogonal geometry option, extension of the 
geometrical capability to treat fully three-dimensional wings (with taper and 
sweep), and allow for transformations to account for boundary layer growth. 
Further, a method was developed to describe the wing surface, and to distribute 
the grid points throughout the domain to allow for the accurate solution of the 
governing equations . In the following section the governing equations and 
computational procedure are described. Thereafter a discussion is presented of 
the new methodology and the calculations that were conducted. 

Computational Procedure 

In describing the overall computational procedure, consideration is given to 
the governing equations: the turbulence model and the numerical algorithm. These 

topics are now briefly discussed. 

Governing Equations 


In the following, the governing equations are nondimensionalized as follows, 

x 1 - with respect to the characteristic length L, the velocity with respect to U*, , 

2 2 

density, pressure and temperature with respect to />a>U and U^, /Cp, 

respectively and time with respect to L/U*, . The viscosity is nondimensionalized 

with respect to . 
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Continuity Equation 


— + — [j/?U k ] = 

31 J [ J ,k 

where J is the Jacobian, p the density, and 
component . 

Momentum Equations 

The momentum equation in the direction is 



where ' ,k' denotes a partial derivative, '|k' denotes a covariant derivative and 

glk a component of the metric tensor . 

In Ref. 12 it was pointed out that the QR Operator scheme requires that the 
governing equations be in quasi -linear form and that the spatial operator in a 
given direction operate on only one variable. For the momentum equation this 
requirement prevents the implicit treatment of certain diffusion terms that arise 
due to the curvature effects. In the usual boundary layer approximations these 
explicitly treated terms would not appear in the equation since they are of order 
0 (Re^/2) or smaller, and should, therefore be of little consequence. 

Since mixed partial derivatives are commonly treated explicitly in orthogonal 
coordinate systems, this same approach is used in generalized nonorthogonal 
coordinates and this concept is extended to include mixed second covariant 
derivatives. All other second covariant derivatives are retained as implicit. 
Since the pressure is specified and impressed upon the viscous layer, its 
specification replaces the normal momentum equation. Thus, the streamwise and 
spanwise momentum equations are the only two retained. A more detailed discussion 
of the derivation of these equations is given in Ref. 12. 


( 1 ) 


u^ is the k ^ 1 contravariant velocity 
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Energy Equation 


For the energy equation, constant stagnation temperature is assumed. 
Neglecting the square of the normal velocity with respect to the squares of the 
other velocity components 


To 


i r 2 x 

T + 2 l U P + W Pj 


912 

h l h 2 


U P W P 


(3) 


where Up and w p are the physical velocity components. These assumptions are 
employed here only for simplification purposes. If warranted, they can be removed 
and the full energy equation can be considered. 


Equation of State 

The equation of state assumes a perfect gas and is given by 



Linearizations 

The following analyses assume a set of linear partial differential equations. 
However, the convective part of the momentum equation and the continuity equation 
are nonlinear, containing terms that involve the product of density and velocity 
components. In order to overcome this difficulty, the procedure described in Ref. 
12 is employed to linearize the aforementioned terms by Taylor series expansion 
about the known time level solution. 

It is important to note that in the governing equations the contravariant 
velocity components are used. However, as noted in Ref. 12, it is advantageous to 
solve for the physical velocity components. Therefore, when the governing 
equations are subsequently cast into a form amenable to the application of the LBI 
scheme, they are transformed so that the physical velocity components appear. 
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Turbulence Model 


In turbulent flow cases, the three-dimensional ensemble -averaged turbulent 
flow equations are considered. The approach taken here assumes an isotropic 
turbulent viscosity, relating the Reynolds' stress tensor to mean flow 

gradients. Using Favre averaging (Ref. 15) the governing equations then are 
identical to the laminar equations with velocity and density being taken as mean 
variables and viscosity being taken as the sum of the molecular viscosity, p, and 
the turbulent viscosity, Pj. 

At this point additional closure assumptions for the Reynolds stresses are 
required, i.e., the evaluation of fi T . There are a variety of approaches 
available, from the simpler mixing length models to the more complicated one and 
two-equation models. Since the method is being applied to wall bounded cases, the 
mixing length model which has worked well in the past for similar flow 
environments (Ref. 16) was chosen. The extension to more complex models could be 
undertaken at a later time if warranted. At that time, the LBI procedure that is 
used for the solution of the momentum equation could be applied to the k and e 

equations . 

Employing the Prandtl mixing length concept, the turbulent viscosity is given 
as 

/i T = pi 2 j* (5) 

where i is the mixing length and $ is the dissipation function, which in 
generalized tensor notation is given by 

$ = | eiieij < 6 > 

As in the Cartesian formulation, * does not automatically reduce to the dominant 
term for standard boundary layers, i.e., 3u/3y in two dimensions and [(3u/3y) 2 + 
(dw/dy) 2 ] 1 / 2 in three dimensions. Hence, provisions are made in the computer code 
that on option retain only the dominant components of the strain which would 
conserve computer time . 

The mixing length formulation is based on McDonald's model (Ref. 17), and is 
given by 
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Z = i^tanh [^-1 D 

oo J 


(?) 


where Z m is the outer layer length scale, and 

D = i — exp(— y + /A + ) ( 8 ) 

where y + takes on its usual meaning. The constants appearing in Eqs. 7 and 8, k , 
\ and A + are .4, .09 and 26.0, respectively, and S is the local boundary layer 
thickness defined as .995 U e . Note that in the limit as y-0, Eq. 7 reduces to 

£ i = kyD 

while for y, large Eq. 7 reduces to 

Z Q Z 

the standard two layer values . 

Numerical Procedure 

The numerical procedure used to solve the governing equations is a 
consistently split linearized block implicit (LBI) scheme originally developed by 
Briley and McDonald (Ref. 15). The procedure is discussed in detail in Ref. 7. 
The method can be briefly outlined as follows: the governing equations are 
replaced by an implicit time difference approximation, optionally a backward 
difference or Crank-Nicolson scheme. Terms involving nonlinearities at the 
implicit time level are linearized by Taylor expansion in time about the solution 
at the known time level, and spatial difference approximations are introduced. 

The result is a system of multi -dimensional coupled (but linear) difference 
equations for the dependent variables at the unknown or implicit time level. To 
solve these difference equations, the Douglas-Gunn (Ref. 11) procedure for 
generating alternating- direction implicit (ADI) schemes as perturbations of 
fundamental implicit difference schemes is introduced in its natural extension to 
systems of partial differential equations. This technique leads to systems of 
coupled linear difference equations having narrow block-banded matrix structures 
which can be solved efficiently by standard block- elimination methods. 
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The method centers around the use of a formal linearization technique adapted 
for the integration of initial-value problems* The linearization technique, which 
requires an implicit solution procedure, permits the solution of coupled nonlinear 
equations in one space dimension (to the requisite degree of accuracy) by a 
one -step noniterative scheme. Since no iteration is required to compute the 
solution for a single time step, and since only moderate effort is required for 
solution of the implicit difference equations, the method is computationally 
efficient; this efficiency is retained for multi -dimensional problems by using 
what might be termed block ADI techniques. The method is also economical In terms 
of computer storage, in its present form requiring only two time -levels of storage 
for each dependent variable. Furthermore, the block ADI technique reduces multi- 
dimensional problems to sequences of calculations which are one -dimensional in the 
sense that easily- solved narrow block-banded matrices associated with 
one -dimensional rows of grid points are produced. A more detailed discussion of 
the solution procedure is discussed by Weinberg and McDonald (Ref. 12) and is 
given in Appendix B. In Appendix A the QR operator scheme is described which is 
used to obtain the spatial approximations. Further details can be found in Ref. 

8 . 


Tasks Considered 

The tasks considered in this contract are now described. 

Geometric Modifications 

This task pertained to the extension of the geometric capability of the 
computer code. In the existing computer code the metric tensor had a special form 
due to the fact that one coordinate is normal to the surface of the body and 
independent of the surface coordinates. Thus, the components g ^3 = g 3 £» for i * 3 
are identically zero. This assumption restricted the ability to efficiently 
resolve boundary layers which have a significant variation of boundary layer 
thickness over the region of interest. If one wished to efficiently redistribute 
grid points normal to the wing by incorporating coordinate transformations that 
account for boundary layer growth, i.e. a y/S transformation, where 5, the 
boundary layer thickness, is in general a function of the surface coordinates, 
then the metric tensor would become full. Since the new normal coordinate will be 
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a function of the streammise and spanwise directions, g 13 and g 23 will be nonzero. 
Furthermore , the variation of the normal coordinate along the surface of the body 
will lead to corresponding nonzero Christo fel symbols. 

To achieve this enhanced generality and efficiency of the computer code the 
geometry arrays were expanded to include the additional nonzero entries and to 
accommodate the three-dimensional variation of the geometric coefficients. With 
regard to the computer code, several terms in the governing equations which were 
excluded previously were added to account for the new geometry. 

Several computations were performed to validate the extended geometric 
capability. At this stage, the two-dimensional case considered in the Phase I 
effort, i.e., Karlsson's experimental data were recomputed using a y/S 
transformation. The y/S calculation considered is again the flat plate turbulent 
boundary layer, but now the outer boundary is permitted to grow at a prescribed 
rate. This allows for greater resolution near the upstream boundary by packing 
more points within the viscous layer. 


Wake Wing Calculations 

In the previous version of the computer code only one surface was treated, and 
that surface was a plane. For real wings, the surfaces are curved and both upper 
and lower surfaces must be treated as well as the wake. Hence, efforts were 
undertaken to modify the code and allow for the consideration of realistic 

geometries . 

The first step undertaken was to define the airfoil shape. This involved the 
specification of a wing cross-section surface, i.e. thickness versus chord or 
distance from the leading edge. Afterwards the arc length was computed and grid 
points were distributed along the surface as desired from accuracy considerations. 
Several different types of airfoils are allowed, including the NACA OOXX series 
and the ONERA type sections. The code is sufficiently general to permit others, 
as well, for which either a formula is prescribed or y vs. x data is given. This 
procedure is employed for the upper and lower surfaces. 

Of comparable importance is the viscous flow in the wake. The wake is 
obtained as an extension of the wing surface. That is, the centerline of the wake 
begins at the trailing edge and extends downstream several chord lengths, which in 
the cases considered, is somewhat further than three chord lengths. The outer 
edge of the computational domain is also prescribed at a fixed distance above the 
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surface, following the shape of the airfoil and for the wake region staying 
parallel to the centerline. 

Two special features of the procedure are of note. First, the trailing edge 
region must be handled carefully. Since the geometric coefficients required by 
the calculation contain derivatives of the metric tensor, the metric tensor must 
remain smooth everywhere in the computational domain. For the trailing edge 
region this means that sharp angles are not permitted. Hence, the trailing 
edge was smoothed and a cusped region (zero slope) was added. Furthermore, the 
local radius of curvature was chosen such that the normal lines eminating from the 
concave surface would not intersect within the computational domain. Second, with 
regard to the wake centerline, in reality it is a double line consisting of the 
extensions of upper and lower surfaces. Special attention was given as is 
described subsequently. 

Figure 1 shows a typical wing/wake coordinate system. The streamwise velocity 
is considered positive in the downstream direction from the leading edge for both 
the upper and lower surfaces. The normal velocity is considered positive directed 
away from the wall. Hence, the normal velocity on the upper surface is positive 
pointing up, while on the lower surface it is positive pointing down. 

The solution procedure is described now. As noted previously, an Alternate 
Direction Linearized Block Implicit method is employed in the solution of the 
equation. In the first sweep, in the streamwise direction, the streamwise 
momentum equation is solved. For three-dimensional flow the spanwise momentum 
equation is also solved concurrently. First, the equations along these streamwise 
coordinate lines are solved for the upper surface and thereafter for the lower 
surface. The streamwise lines in these two regions extend from the leading 
boundary to the outflow or downstream boundary at the termination of the wake. 

This leaves only the wake "double line", which must be solved. The line extends 
from the trailing edge to the outflow boundary. Only one of these lines is 
computed, and the values are set on the other line. This completes the first 

sweep . 

In the second sweep the equations are solved on the normal lines to the 
surface (cf. Figure 2). First the top surface and then the lower surface regions 
are computed until the onset of the wake region. In the wake additional 
manipulation of the equations are now performed since the entire region as a whole 
must be solved from the bottom to the top. This manipulation involves ordering 
the equations appropriately and accounting for the double line on the centerline 
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of the wake. At the boundary layer edges, above and below the airfoil, the 
streamwise velocity is specified as the boundary condition. On the wing 
surface no slip is specified, i.e. zero velocity. For three-dimensional flows, 
the spanwise edge velocity conditions are also specified similarly to the 
streamwise conditions at the boundary layer edge. 

Different forms of boundary conditions are specified for the continuity 
equation. Since the continuity equation is discretized as a two point trapezoidal 
integral form, no boundary conditions as such are specified at the boundary layer 
edge but rather, the governing equation is solved there. However, at the 
centerline an additional condition is required. Since the velocity is not known 
there, a priori , it obviously cannot be specified. Instead, a smoothness 
condition is enforced, i.e., the second derivative of the normal velocity is set 
to zero. Note that for the momentum equation the governing equation itself is 
solved at the centerline. 

Once the two sweeps are completed the velocity components are transformed into 
their physical normal to the wall and streamwise components. 

Two-Dimensional Steady Turbulent Flow for NACA 0012 Airfoil 

The case performed is the turbulent flow over a NACA 0012, symmetric airfoil 
at a 4.86° angle of attack. The frees tream chord Reynolds number was .48E+07 
and the mean freestream Mach number 0.599. The freestream temperature was assumed 
to be 300° K. The computational domain was chosen with inflow boundary located at 
x = .1 ft (.1 chord) and the outflow boundary located at x = 4.6 ft, while the 
outer edge was set to a constant value of .25 ft in the direction normal to 
airfoil surface for all x (the airfoil leading edge corresponds to x = 0 and 
trailing edge corresponds to x = 1) . The grid distribution in the normal 
direction is based upon a hyperbolic tangent function. In the streamwise 
direction, a distribution based on series of error functions was used (Ref. 14). 
There were 49 mesh points distributed in the normal direction and 118 mesh 
distributed In the streamwise direction around the airfoil (62 grid points in the 
body region and 56 grid point in the wake region) . From the pressure 
coefficients, the freestream velocity was calculated using isentropic 
relationships . 

At x = .1, the inflow boundary, the displacement thickness of the upper and 
lower surfaces are assumed to be the same as that obtained from the integral 
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method, which are .60207E-03 ft and . 33634E-03 ft, respectively. 

The displacement thicknesses Reynolds numbers, Re£ , at the above locations 
are therefore, 2800 and 2000, respectively, which can then be used to calculate 
the corresponding Cf via the Clauser formula 


J 2/C f = 5.6 log R eS* + 4.3 

The Cole's law velocity profile, i.e.: 


u 1 1 
— = — In 

r yu T i 

1 „ 27t(X) . o 

| + C + — - sin^ 

1 X 

u r 

T K 

V 

i K 

2 ^ 


Wall -Wake Law 


Viscous Sublayer 


(9) 


where 




U T = J T w /p 


and tt(x)//c is evaluated from the condition that u = at y/8 = 1; furthermore, 
constants k and C are set at .41 and 5.0, respectively. 

As described in Ref. 15, the relationship among ReS*, 8 , Cf, II (x) is 


k (R e5 * -65) 

SV T /v 


i + n (x) 


Eliminating II (x) from Eq. (9) and Eq. (10) yields 


( 10 ) 



- In ) + ^ \-ReS* - 65 ]- - 
K X l J K 


(ii) 


where 


X = 8\X r /v 

The 8 value, therefore, can be obtained from Eq. (11) by a Newton Ralphson 
iteration procedure which, together with Cf determines the required velocity 
profile. Once the streamwise velocity profiles are obtained the normal velocity 
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can be determined by assuming constant density for the purpose of setting upstream 
normal components and then integrating the continuity equation. 

In the streamwise direction, the boundary layer option was employed, that is 
the streamwise diffusion terms were ignored and a backward difference 
approximation was used for the streamwise convective terms. The boundary 
conditions stipulated on the body surface were no-slip and zero normal velocity. 

At the outer edge of the viscous layer the magnitude of the streamwise velocity 
component was also prescribed. The value of the normal velocity component is not 
set, but rather computed as part of the numerical solution, as is the practice in 
standard boundary layer procedures. At the inflow boundaries, velocity profiles 
are fixed. Furthermore, the intermediate boundary conditions employed on the 
first sweep are the physical ones. For steady problems, the imposition of 
physical intermediate boundary conditions did not impair the quality of the 
solutions obtained. These results are in keeping with the analysis of McDonald 
and Briley (Ref. 15) for second order spatial scheme. A comparison between 
results obtained from VISTA 3-D and an integral method provided by NASA LARC 
personnel is shown in Figs. 3 and 4. Fig. 3 shows the displacement thickness as a 
function of distance along the airfoil. As can be seen, the agreement of 
predicted value between VISTA 3-D and the integral method for pressure side of the 
airfoil (lower surface) is excellent. For the suction side of the airfoil (upper 
surface), these two results were not in as close agreement, especially when x is 
large . 

Similar curves for momentum thickness are shown in Fig. 4. Again, for the 
lower surface, the agreement is excellent; however, for the upper surface there 
was some disagreement between the two solution procedures, with the integral 
method overpredicting the values relative to the VISTA3D code. 

Two-Dimensional Unsteady Turbulent Flow for NACA 0012 Airfoil 

For the unsteady calculation the sample problem provided by NASA LARC 
personnel was a NACA 0012 symmetric airfoil which performs a sinusoidal pitching 

o 

oscillation at the frequency of 4.789 Hz with amplitude of 1 . The frees tream 
chord Reynolds number was .48E+07 and the mean freestream Mach number was .599. 

The procedure for obtaining unsteady flows is similar to the steady solution 
procedure. In the steady state case, the inflow and outer edge boundary 
conditions are prescribed to be invariant in time. For the unsteady case. 
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however, the velocities at these boundaries are allowed to change in time. At any 
instance the outer edge velocity is determined via the given pressure coefficient 
obtained from an inviscid calculation. In specifying the upstream velocity 
profile the present procedure adopts the same approach as does the steady state 
case. That is, the upstream profile is still given by the Cole's wall-wake 
velocity profile, but now C f and S* are allowed to be periodic functions of time; 

i.e. , 

C f = c f © ( 1 + ACf cos(wt + * Cf )) ( 12 ) 

where Cf and S* are the mean (time averaged) skin function and displacement 

thickness, Cf^ and 6* are their respective amplitudes of oscillation and <*> ^ and 

<p gie their respective phase shifts. This procedure introduces additional unknowns. 

For the mean quantities C f and S* , the steady state values are used. While the 

o o 

other four quantities are determined by the characteristics of the current 
problem. For this problem the oscillatory frequency, w , is 4.789 Hz, the chord is 
1 ft and the freestream velocity is 694 ft/s. Hence, the corresponding reduced 
frequency, based on chord, is .0433. This low reduced frequency implies that the 
quasi-steady flow is a valid assumption; physically this means that oscillatory 
changes are much slower than convection changes. For this reason, the phase angle 
<f> cf and UT<p s * in Eq. (9) would be insignificantly small and therefore is set to 
zero in the current study. Due to the quasi-steady characteristic, given the edge 
velocity in Cole's velocity profile the instantaneous C f and S would then be 
determined via the Newton-Ralphson iteration. In this study, however, a reasonable 
oscillatory amplitude would be given to Cf and S* to investigate the unsteady 
phenomena . 

However, the low reduced frequency nature of the unsteady problem presents 
some constraints for the computations. First, the temporal discretizations are 
determined by the smallest time scale of the problem, i.e., the convective 
characteristics of the flow field. In general this means that for each time step 
that a particle in the freestream should not travel more than 10%-20X of the chord 
length. Based upon this criterion, approximately 1500 time steps are required in 
order to achieve temporal accuracy, significantly more time steps than were needed 
in the calculation of Ref. 6. As noted previously, 118 points were used in the 
streamwise direction. In order not to use an extensive amount of points in the 
trailing edge region, the inviscid pressure distribution was slightly modified to 
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diminish the large gradient there. This permitted the use of a coarser grid. 
Physically, one expects the inviscid pressure distribution to be higher than the 
viscous distribution, so that the aforementioned modification is not unreasonable. 
Thereafter, several test runs were performed to determine the optimum time 
increment that would not compromise the computational efficiency and numerical 
stability. The final selection of At is .00029 sec, which corresponds to 
nondimens ional time steps of .208 (i.e. for each time increment a freestream 
particle would travel about 20X of the chord length) and 720 time steps per cycle. 

It is noteworthy that in the Karlsson calculations (Ref. 6), due to the low 
freestream Mach number, the reduced frequency based on chord length was in the 
order of 1 to 10. For this reason, the time step could be determined from the 
external oscillatory frequency and therefore, 36 steps per cycle in most cases 
would produce satisfactory results. For a typical helicopter blade, the reduced 
frequency under normal operating condition is usually very low, hence the above- 
mentioned computational concerns would also be encountered. For comparison 
purposes a turbine blade, has a reduced frequency in the order of 1, and thus the 
aforementioned problems would not appear at all. 

In presenting the results, the skin friction, C f , the displacement thickness, 
g* and the streamwise velocity profile at the measuring station were Fourier 
decomposed Into their harmonic components. 

f (r ) = ^ + l ja n coswnt + b n wnt| 


where 


tj+T 

a Q = - f f (£ ) coswn£d£ 
it . 


t-j+T 

a n = - J f(£)coswn£d£ 

7T 

t i 


tj+T 

b n = - J f (£ ) sintt>n£d£ 

It 

^1 
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and where t^ is the time at the start of the integration, and the period T = ln/u). 

Although the code allows for the determination of any number of harmonics, 
only the first two were obtained. Therefore, the Fourier series representation 
mean velocity is a Q /2 , the in-phase component of the first harmonic is a^ and the 
out-of -phase component of the first harmonic is -b^. For the evaluation of the 
Fourier coefficients Simpson's integration scheme was used and all data points for 
the second cycle were sampled. 

Table 1 is the result of the skin friction coefficient. It is seen from this 
table that the symmetry conditions (in the periodic sense) of the upper and lower 
airfoil is maintained, which is a necessary condition for the validity of the 
results. It is also seen that the mean skin friction coefficient is about fifty 
to one hundred times larger than the first harmonic oscillatory amplitude and the 
first harmonic amplitude is an order of magnitude larger than the second harmonic 
amplitude, which implies that for this particular problem the zero oscillatory 
amplitude of the friction coefficient for the upstream velocity profile is a valid 
assumption. Another interesting point is that the phase angle changes along the 
airfoil stations, which indicates that care must be taken to determine the 
upstream profile for a more general problem. 

Table 2 and Table 3 present the Fourier analysis of the displacement thickness 
and momentum thickness, respectively. The general observations are similar to 
those of the skin friction coefficient, i.e.: 1) the symmetry conditions are 

satisfied; 2) the mean value is much larger than the oscillatory amplitude and the 
second harmonic oscillation is negligible; and 3) the phase angle is a function of 
location along the airfoil. 

Table 4 shows the Fourier analysis of the streamwise velocity component at 
various stations along the normal direction at x = .8 (chord). It can be seen 
from this table that the mean value is the most important Fourier coefficient and 
therefore there is a significant phase shift from the wall to the outer edge. 

The Fourier analysis of the pressure coefficient obtained from an inviscid 
analysis is shown in Table 5. As can be seen in this table the mean value and 
first harmonic amplitude of the pressure coefficient at the corresponding stations 
on the upper and lower surfaces are practically the same and the oscillatory 
angles are nearly 180 degrees out of phase. Furthermore, the amplitude of the 
first harmonic oscillation is about fifteen to twenty times larger than the second 
harmonic oscillation, indicating the first harmonic dominancy of this particular 
problem due to the small pitch oscillation. 
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v/<S Calculation 


In order to test the nonorthogonal capability in two-dimensional and 

three-dimensional flows, several demonstration calculations were conducted. For 

the two-dimensional case, a y/S calculation was undertaken for the incompressible 

turbulent flat plate flow. In previous calculations a cartesian coordinate system 

was employed, consisting of streamline and wall-normal coordinates. Since the 

boundary layer grows rather rapidly from leading to trailing edge; near the 

upstream boundary there would be very few points contained within the boundary 

layer. In order to cluster more points into the boundary layer, the outer edge 

u 

was set at twice the local boundary layer thickness and allowed to grow at a 
rate. Thus, the 'streamwise' coordinate lines were no longer parallel to the wall 
but were curved, while the normal to the wall coordinates remained unchanged. 
Hence, the intersection of these two coordinates led to a full nonorthogonal 
coordinate system in which all geometries terms would be tested. It should be 
noted that in order to apply the method to generalized coordinates not only must 
one obtain the appropriate geometric coefficients, but one must compute the 
applicable velocity components. Since the governing equations are solved along 
coordinate lines, the physical velocity components must be transformed correctly 
to account for the curvature of the lines. Further, the pressure gradient term 
which is imposed on the boundary layer, at the outer edge was also accounted for 
to assure that one of the principal assumptions that the pressure remain constant 
through the layer (normal to wall) be enforced. 

The results of the calculations for the nonorthogonal cases were compared to 
the cartesian case, with regard to velocity profile, skin friction coefficient 
displacement thickness and momentum thickness. These comparisons indicated that 
the results were well-behaved and given indistinguishable values, verifying the 
procedure. These modifications were then employed to compute the flow over a 
NACA 0012 airfoil. 

Three Dimensional Unsteady Flow over a Skewed Nonorthogonal Coordinate System 

For the three-dimensional calculation, the other aspect of the nonorthogonal 
coordinate system was investigated. In this case the nonorthogonality was limited 
to the surface, in the streamwise and spanwise directions instead of the 
normal/s trearawise direction considered previously. The case undertaken was the 


- 19 - 



three-dimensional flow over a flat plate in which the streamwise flow is skewed to 
the leading edge at 45°. Fig. 5 shows the orthogonal case, which was studied in 
Ref* 6. In figure 6 is shown the case considered in the present study in which 

o 

the streamwise coordinate lines are skewed at a 45 angle to the spanwise lines, 
and the normal lines to the surface remain normal. In this coordinate system all 
lines remain straight. The resulting spanwise velocity components in the new 
nonorthogonal system should remain zero. This indeed was the case for the 
calculation* The calculation was run in the steady and unsteady modes and the 
results compared well with the cartesian calculation. 

Three-Dimensional Wing Studies 

The first step in treating the flow over actual three-dimensional wings is to 
specify the geometry. This entails, as with the two-dimensional airfoils, the 
description of the wing cross-sectional shape, the specification of inflow and 
outflow boundaries and the boundary layer edge. A similar procedure to that which 
was used for the two-dimensional case was considered here for three -dimens ions 
for specifying the surface and distributing the grid points. Please refer to that 
section for details. The one important difference for the three-dimensional case 
is that now, in addition to the streamwise and normal coordinate distribution, a 
spanwise distribution is required. Such a distribution is obtained by considering 
the wing consisting of a sequence of spanwise 2-D sections. The number is chosen 
from a consideration of the variation of spanwise shape; i.e., taper and sweep and 
the requirement of smooth geometric coefficients. Hence, the three-dimensional 
flow field is built up from a series of two-dimensional sections. An example of 
such a coordinate system is given in figures 7 to 13, in which an ONERA wing is 
shown. 

Figures 7 and 8 are 2-D sections of the wing. 

Figures 9, 10 and 11 are 3-D perspectives with spanwise coordinates on the 
wing and streamwise normal coordinates at the tip and root sections. 

Figures 12 and 13 show the spanwise coordinates on the wing's surface. Note 
that the wing terminates at a fixed streamwise location. 

Note that at the coordinate system begins .1 chord downstream of the leading 
edge where the inflow boundary conditions are specified. 

The solution procedure is similar to the two-dimensional calculation with the 
exception that the spanwise momentum equation must be solved. This leads to a 
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three sweep method (see Appendix) with an increase in matrix block size. 

As can be realized, the inclusion of a third dimension increases the 
complexity of the geometric specification since first and second variations of the 
metric need to be considered in three dimensions. For the skewed boundary layer 
described previously, although the streamwise and spanwise coordinate system was 
nonorthogonal , there was no curvature. This reduced the number of nonzero 
geometric terms. However, in the current cases where there is curvature of the 
coordinate lines , additional terms now present themselves . The entire geometry 
array was thus carefully checked for accuracy and smoothness . In order to verify 
the solution procedure a test case was constructed in which the entire 
three-dimensional flow code would be exercised, but in which the flow remains 
essentially two-dimensional. In this case the wing consisted of 5 identical 
spanwise planes (2-D airfoil sections), as shown in figures 10-13. First, the 
two-dimensional counterpart was solved in the steady state. This solution was 
used as the inflow boundary condition. At the other spanwise boundary (plane 5) 
an outflow boundary condition was set so that the first derivative velocity 
components was equal to zero. This boundary condition transforms the flow field 
to one which does not vary in the spanwise direction, thereby retaining a 
two-dimensional character. Calculations were conducted for this test case. 

Although the flow field appeared to be two-dimensional through most of the 
flow field, there were regions where anomalous velocities appeared. Efforts were 
undertaken to discover the source of this problem. Unfortunately these efforts 
were not totally successful. Areas in the code were identified which could 
contribute to the observed results. In particular the computation of the source 
terms which arise from the evaluation of terms at the lagged time step were 
identified as the prime source. However, further work could not be conducted 
under the present contract. The results obtained to date are very encouraging and 
indicate that calculations of the type considered are realistic and attainable. 


CONCLUSIONS 

In this report, a method for solving unsteady flows over two- and 
three-dimensional wings was described. The initial computer code was extended to 
treat three-dimensional geometries in a nonorthogonal coordinate system containing 
coordinate lines with curvature terms. The entire wing is solved as a whole, 
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including the upper and lower surfaces and wake. The method was tested in two 
dimensions for the unsteady flow over a NACA 0012 airfoil. In three dimensions, 
the coordinate system, geometry and governing equations were extended to treat 
realistic wings. Although complete three-dimensional solutions were not obtained, 
this effort has laid the groundwork for the computation of such flow fields of 
interest . 
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appendix a 

SPATIAL DIFFERENCE APPROXIMATIONS 


QR Operator Notation 

In this section, implicit tridiagonal finite difference approximations 
to the first and second derivatives and to the spatial differential operator 
are considered. The QR operator procedure for generating a variety of 
spatial discretizations is also introduced. As special cases,, standard 
second-order finite differences, first-order upwind differences, fourth-order 
operator compact implicit (OCI), fourth-order generalized OCI and exponential 
type methods are obtained. Since all these schemes are of the same form 
(cf. below), a single subroutine which defines the difference weights is all 
that is required to identify the method, while leaving the basic structure of 
the program unaltered. The rationale for the use of the QR approach in the 

present problem is discussed in detail in Ref. 8 • 

The QR fore.el.tion .lions for MU methods and pers.it. the treatment of 

systems of coupled equations. i.e., LB1 methods. Although v.rUble mesh 
schemes e.n be employed within the QR fra»««ork, it is believed preferable to 
use analytic transformations to obtain a uniform computational mesh, 
attention is restricted to uniform mesh formulations. 

The general concepts and notation will be introduced for two-point 
boundary value problems and then the methodology will be extended to more 
general linear and nonlinear parabolic partial differential equations in one 
dimension. The application of QR operator method to multidimensional 
problems is discussed in the section pertaining to the LBI scheme. 

Consider the two-point boundary value problem 

Liu) - o( x )u xx + b(x)u x + c(x)u - fix) (b-1) 


with boundary values u(0> and ull) prescribed. Derivative boundary 
conditions, although not discussed here, can easily be incorporated into the 
framework of the QR operator notation. Let the domain be distressed so that 

xj = (j-1 )h , j= 1 , 2 


. , J + 1 , and Uj 


u(xj) , Fj 


^xx ( x j ) and h = 1/J is Che mesh widch. The numbering 


convention was chosen here Co be compatible with FORTRAN coding. 
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Without loos in generality for .<*> * 0. Eq. (A-l) can be divided by 
,(,) to that we may treat instead the following equation 


L(u ) • u„ + b(x)u x + c(x)u - fix) 


(A-2) 


where 


blx) - b(x)/?(x). c(x) - ZM/oM ond fix) - f (x)/a(x) 


Subatituting the finite difference approximations to the first end 
second derivatives 

D "" ^i— I # . \ . 


'O , Jill 
r 2TT U J “ 2h 


Ui+ « " U hL . F . + 0(h‘) 


DD Ui_,-2Uj+U: + i \.rhfK 2 ^ 

-^U, ~ ' J — ' S |- U xx (,t i , + 0(h 1 


(A-3) 

(A-A) 


h 2 

into Eq. (A-2) and rearranging, we obtain j 


L(u) 


[-^r " lEir] °H + h " tN° j + [ "S* - + Zh - ] ° J *' ‘ '* 


( 1 A) 


or 


['- ^]°H + l h V 2 l u l + [ , + _ i L ] U i 


(A-5) 


)♦' • h2 'i 


where RCj - hbj is the cell Reynolds number. 

Equation (A-5) can be generalised by introducing operator format, i.e. 


rjUj., d-r'uj + rfu,,, - ht <dj fj-s + q j *, + 


(A-6) 
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where the superscripts (— ) roinu6, (c) center, and (+) plus indicate the 
difference weight that multiplies the variable evaluated at the (j-1), (j) 
and (J+l) grid points respectively, and where the rj's and qj's for grid 
point j are functions of h, bj_j, b j » bj + j, c j-l* cj and Cj + 1* 

Comparing Eqs. (A-5) and A-6) we can identify the r^’s and q^’s, viz.. 


rj - I - Rcj /2 qj “ 0 

r j - h 2 Cj - 2 Qj “ I 

r J - I + Rcj /2 qj " O 

We now define the tridiagonal difference operators Q and R 

«[«,] - r j u i-t +r i u i + r I u J*> 

°[ f j] ' q l f |-' + q ) f i + q l f l« 


(A-7) 


(A-8) 


Noting that L(u) - f and substituting Eq. (A-8) into Eq. (A-6), we obtain 

r[ u j] - h 2 o[f J ]« h 2 o[ L(u)j] 


Alternatively by employing the inverse operator Q 1 an expression for 
L(u)j can be obtained 


L(u) j 



CT'RUj 


(A-10) 


- 28 - 



For standard central finite differences Q = Q 1 = I, the identity 
matrix, so the spatial operator can be given explicitly in terms of Oj_j, 

Uj and Uj + |. However, in general, for higher order methods whereas Q is 
tridiagonai Q -1 is a full matrix, and the spatial operator cannot be given 
explicitly in terms of the variables at adjacent grid points. Hence, Eq. 
(A-10) provides a method for expressing the spatial operator for a wider class 
of difference approximations. The formalism in Eq. (A— 10) is also applicable 
for first and second derivatives appearing alone Ccf. Ref . _18) . In Refs. 8 
•arid 19 a technique *lue to Berger, et al/'ls described far constructing fourth 
order tridiagonal methods which possess a monotonicity property as the cell 
Reynolds number is increased, R^ ♦ **. This type of scheme is an option in 
the computer code. 


APPENDIX B 

LINEARIZED BLOCK IMPLICIT SCHEME 
Consider a system of nonlinear partial differential equations 

a 5. = <b-« 

' 1 

where t is a vector of unknowns and I is a source term vector which is a 
function of x 1 , x 2 , x 3 and t. Extension to source terms which are functions 
of I are discussed in Ref. 15 . is a three-dimensional nonlinear 

differential operator and the matrix A appearing in the momentum equations is 
equal to pi where p is the density and I the unity matrix. 

Equation <B-1) may be centered about the n+B time level, i.e. t n+B = 

(n+B )At = nAt+BAt = t n +&At, and written 

a -/3 [ /A, v"' 0 < B - 2 > 

where 0 < 0 < 1 is a parameter allowing one to center the time step, i.e., 

0=0 corresponds to a forward difference, 0 = 1/2 to Crank-N icol son and 0 = 

1 to a backard difference. 
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After Linearizing Eq. (B-2) by Taylor series expansion in time about the 
n £ h time level by the procedure described in Ref. 15 to give a second-order 
linearizat ion , we obtain 

A "[ &"]/&, , y n [* n+ ^ n ] - 3 > n &" + <»-’> 

where X is the linearized differential operator obtained from 

The difference between the nonlinear operator and the linear 
operator X is .defined as bP 1 ** — X n * At the intermediate level 

n + 6, $n+0 is represented as 

= /? <f> nfI + (!-/?)$>" (B “ A) 

Using these relationships and dropping the vector superbar for convenience a 
two- level hybrid implicit - explicit scheme is obtained 

A"(4» nt, -4. n )/ai = ^ n (4. n+, -4>") * M n <t> n + y n+e ( B -« 

The vector represents all of the terms in the system of 

equations which are treated explicitly. More about this will be said later, 
but for the moment note that may be approximated to the requisite 

order of accuracy by some multilevel linear explicit relationship, or 
approximated by tj* n with a consequent order reduction in temporal accuracy. 

The operator X is now expressed as a sura of convenient, easily 
invertible suboperators X ** X \ + X 2 + • • * • X m » In the usual ADI 
framework these suboperators are associated with a specific coordinate 
direction. Further, it is supposed that these suboperators can be expressed 
in the QR notation introduced earlier- Writing and M n 4> n as a 

single source term S n ‘ f ^, Eq. (B-5) is written as 

/ A! [<J>"* L <J> n ] + S n+/3 (b-6) 
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. . col it into a sequence of easily 

t i system efficiently it is 6 P lxt 

° 60 VC 1 {ration of the procedure of Douglas 

invert ible operation, follo-log a gener.lrr.t.on f 

an, Cunn (Ret. 18) in i» n.tur.l exten.ion to systems o£ p rt 1 
differential e,natio„e. The Douglas-*™ aplttttng of Eg. (B-6) 
the following three-step procedure 

A n[ *-.*-]/£, + [A n+ A" + A n K +s 

A “[*“ «T)/a, = - &?[*“- *1 + + s ** 

A n [$*”$"]/fit = &:[*' -*"] +/34 n [«>*' -< !>n ] + pt'A® - <p 1 

+ \j” + s" +f> «- 7 > 

which can be transformed to the alternative form 


(B-8) 


[a- -A./ wr] [**-*"] • Ai[./>^ n K+A's n+/J 

[a"-A./ 34"] = A'[<t»~-<f>"] 

= A"[4>**-4> n ] 


If the intermediate levels are eliminated, the scheme can be written in the 
so-called factored form 

(A"-/3At^ 1 n XA n r(A"-)9Ati , 2 n XA n ) (a°-/?A1-4 ><* -<!>)= <B 9) 

8 K/;tA"fA n > 5 ' %4,s ^ 

The Obi formulation given in Eg. (B-8) Is directly applicable for A 

„ a in 0 _1 R operator format. Consideration of 
operators represented in Q P ror 

p . .... nnri the removal of the inverse operator 

intermediate boundary conditions and the remo 

Q-l is given in Ref. 12. ^ / can be spl i c into any 

lt is worth noting that the operator 3> or X P 

„ Hirh need not be associated with a particular 
number of components w . tbe 

. - r • _ As pointed out by Douglas and Gunn (Ref. 11), 

coordinate direction. A P<-> . • „ hv 

, i «- that the associated matrices oy 

criterion for identifying sub-operators is that the 

„ , pc ," (i e narrow-banded). Thus, mixed derivatives and 

easily solved ^l.e., 
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complicating terras which might inhibit the use of OCI can be treated 
implicitly within such a framework, although this would increase the number 
of intermediate steps and thereby complicate the solution procedure. 

An inspection of Eq. (B- 8 ) reveals that only the linearized operators 

• j£ n *nd P n appear. Indeed, the computer code employs this feature by 
12 3 

evaluating these three operators before the first sweep, storing them and 
accessing them as needed in the subsequent three sweeps. In addition, the 
terms arising from the nonlinear terms are immediately absorbed into S n+ ^ 
as they appear, allowing for an efficient evaluation of the terms, in the 
differential equations. 


The spatial operators appearing in the differential equations 

X j, an< * -^3 must be identified at least formally in order to isolate 

the coefficients that are to be used in the construction of the Q and R 
operators. These operators can be represented in standard form at each grid 
point, i.e.. 




(B-10) 


In Eq. (B-10) the first subscript of $ indicates the velocity component 
(associated with the corresponding direction) and 11 , 11 indicates a 

derivative. The subscripts of the a?^ refer to the direction (i) and the 
term in the equation (j) respectively. Note that the equation is in 
quasi-linear form, since the coefficients of the derivative operators need to 
be identified, for use with the QR operator technique employed here. 

Alternate schemes have been proposed by Leventhal (Ref. 20) for equations in 
conservation form but are not considered here. In the following section, a 
description will be given of how this entire operator is discretized by 
employing the QR operator format, and how the discret izat ion is incorporated 
into the LBI framework in order to solve the system of equations (B-8). 

The continuity equation is considered first. Since it is a first-order 
partial differential equation it does not have the standard form of Eq . 

(B-9) . Furthermore, in the linearization process p has been eliminated in 
favor of the u 1 velocity components so that the continuity equation has 
become an equation for the three velocity components, and not density. 
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An inspect ion of the system of emotions under consideration reveals 
that substantial savings can be realized if the equations are partion.d 
appropriately. Due to the use of a boundary layer coordinate system, the 
normal velocity appears only in conjunction with terms associated with the 
normal "3” direction in the two momentum equations. Hence, in the first two 
sweeps where directions “l” and "2" are implicit one is required to solve 
only for the two corresponding velocity component, in the streamwise and 
spanwise momentum equation without the need of considering the continuity 
equation. However, on the third sweep where all 3 velocity components 
appear, one must solve all 3 equations. This strategy reduces the solution 
procedure to the inversion of two 2 a 2 block matrices and one 3 x 3 block 
matrix rather than three 3x3 block matrices which leads to a substantial 
reduction in computation time. If the full Navier-Stokes equations were 
considered (including a normal momentum equation) the aforementioned 
partioning could not be applied since the normal velocity would appear in all 

three sweeps. 

The question that arises is how to appropriately split the continuity 
equation, since it need only be solved on the third sweep. Here again the 
Douglas-Gunn formulation leads to the appropriate choice. The continuity 
equations written in conservation form is. 


dp_ J_ 

at j 


dx 


r [ J / >uI ] = 0 


(B-ll) 


After linearizing 


and eliminating P, the increment form is obtained 


A "Au"*' + B n d/*'+ 22 ^[,VAu“f."B"dx W f; 4 .1 


At d 

j a? 


[ J/)U ']" + [(/ + „"A n )Au"« + ( U V)Aw"-'] 

+ 22. A- r ( " + w°B n ) Aw n+I + (w n A n )Au n+ ^] 

J dx z L r 


(B-12) 


where all the velocity components are the coot ravar iant components u - u , 


w = u 2 and v = v 3 . J is the Jacobian and 
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8 r 



S„ U " + g *2 

g 22 w n + g l2 


n 


W 


n 


U 


] 

] 


By employing the Douglas-Gunn procedure, Eq. (B-12) is represented as a 
third sweep equation, and a consistent approximation is obtained to the 
continuity equation, i.e., the x 1 derivative term is evaluated at the * level 
and the x 2 derivative term is evaluated at the ** level. The values of the 
intermediate derivative terms are obtained after the solution of the first 
two sweeps of the two momentum equations. Note that these terms do not 
contain the normal velocity. The equation can thus be written in symbolic 
form 


_n . n+i . „n A .n+i^ 

A Au + B Aw . + — 


d 

dx 5 


[ J ^A° v°Au 


n+, + v n B n Aw n+, + 


n . n+i 
o Av 


}] 


(B-13) 


= S’ 




a 

*aF 


[*{}]*-!>“■& r-{ >] 


Since the only term involving v is in the x derivative term, one can 

. 3 

directly integrate the equation with respect to x , i.e. 


(A 


A"Au"« 


b°Aw 


"*'1 d *' 



. n . n+1. 
A Au + 


n _n . n+t 
v B Aw 


n A n+i 
o Av 





tfAt.r j** 

tN 


The next section describes how this is done very easily via the QR operator 
scheme. The concept of integrating directly the continuity equation is not 
new. Davis (Ref. 21) in his coupled procedure for Che solution of 
two-dimensional steady boundary layer equations used a trapezoidal rule to 
integrate the continuity equation. Weinberg (Refs. 22 and 23) also used a 
fourth-order Simpson integration scheme to solve the compressible boundary 
layer equations. Such procedures are, stable and offer a viable alternative 

to approximating the derivatives by finite differences. Note Chat 
conceptually the continuity equation in integrated form is treated on each 

sweep of the Douglas-Cunn splitting, although in actuality this can be viewed 
as having the same form as each sweep and the integration operator can be 
incorporated into the _/ and difference operators, and as a result the 

stability and consistency of the original splitting is retained. 
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X Location (ft) 

AO 

A1 

B1 

A2 

B2 



.38967E-02 

.67699 E-07 

-.22287E-06 

.4S814E-08 

.45872E-08 

a> 

o 

Ctf 

t 

■ 

•42797E-02 

-.17854E-05 

.1071 5E*03 

.17208E-05 

.14361 E-05 


.37928E-02 

.18962E-04 

.58182E-04 

.22280E-05 

.28273E-05 

r® V.' ; . 

.33806E-02 

.25652E-04 

.34367E-04 

.16299E-05 

.21 94 IE-05 

3 

(0 

0.46 

.30506E-02 

.31949E-04 

.12734E-04 

.189S3E-05 

.28714E-05 

w» 

0.55 

.27747E-02 

.33389E-04 

-.39809 E-06 

.17S04E-05 

.28084E-05 

<D 

CL 

0.64 

.25333E-02 

.33425E-04 

-.1 1757E-04 

.57673E-06 

.10136E-05 

CL 

z> 

0.73 

.23010E-02 

.31304E-04 

-.20715E-04 

.88873E-06 

.16173E-05 


0.82 

.20370E-02 

.26728E-04 

-.271 19E-04 

.54921 E-06 

.11568 E-05 


0.91 

.16355E-02 

.22326E-04 

-.36889 E-04 

.54578E-06 

.13082E-05 



■SB 

.38967 E-02 

-.65485E-07 

.22272E-06 

.51117E-08 

.20698E-08 

o 


.42784E-02 

.14041 E-05 

-.10725E-03 

.17291 E-05 

.26166E-05 

o 

3 


.37919E-02 

-.19271 E-04 

-.58279E-04 

.23958E-05 

.34449E-05 

l_ 

3 


.33800E-02 

-.26004 E-04 

-.34382E-04 

.18662E-05 

.24906E-05 

<0 

0.46 

.30502E-02 

-.32274E+04 

-.12727 E-04 

.21981 E-05 

.29963E-05 

k_ 

<D 

s 

0.55 

.27745E-02 

-.33953 E-04 

.47986E-06 

.20226 E-05 

.27750E-05 

0.64 

.25332 E-02 

-.33865E-04 

.1 1873E-04 

.8341 4 E-06 

.85366E-06 

3 

0.73 

.2301 0E-02 

-.31758E-04 

•20835E-04 

.11379E-05 

.13654E-05 


0.82 

.20372 E+ 02 

-.27238E-04 

.27251 E-04 

.75024 E+06 

.85661 E-06 


0.91 

.16360E-02 

-.22822E-04 

.37076E-04 

.70514E-06 

.91072E-06 


Table 1. Fourier coefficients of the skin friction coefficient, 

M = .599, a = 0°, Re = .48 X 10, w = 4.789 Hz, Amplitude = 






X Location (ft) 

A0 

A1 

61 

A2 

B2 


0.10 

.37558E-03 

.10251 E-05 

-.34408E-05 

.41487E-07 

.37172E-07 


0.19 

.60572E-03 

.90040E-05 

-.14388E-04 

-.7221 2E-06 

-.771 57E-07 

o 

0.28 

.88322E-03 

.16692E-04 

-.28909 E-04 

-.1 0799E-05 

-.24535E-06 

w 

t: 

0.37 

.1 1733E-02 

.23646E-04 

-.42388 E-04 

-.12803E-05 

-.29465E-06 

13 

CO 

0.46 

.147S3E-02 

.31001E-04 

-.57202E-04 

-.16526E-05 

-.62924E-06 


0.55 

.17918E-02 

.37962E-04 

-.72360E-04 

-.18996E-05 

-.84942E-06 

CD 

Ol 

0.64 

.21279E-02 

.44949E-04 

-.89561 E-04 

-.17706 E-05 

-.53297E-06 

CL 

3 

0.73 

.24990E-02 

.51640E-04 

-.10957E-03 

-.23653E-05 

-.11 151 E-05 


0.82 

.29499E-02 

.58006E-04 

-.13420E-03 

-.27464 E-05 

-.13336E-05 


0.91 

.36501 E-02 

.68505E-04 

-.17760E-03 

-.36737E-05 

-.21 136E-05 



0.10 

.37562E-03 

-.10288 E-05 

.3437 9 E-05 

.48457E-07 

-.26746E-08 

o 

0.19 

.60583E-03 

-.90139E-05 

.14406E-04 

-.66228E-06 

-.24148E-06 

o 

CO 

0.28 

.88345E-03 

-.1671 IE-04 

.28933E-04 

-.96239E-06 

-.57742E-06 


0.37 

.1 1737E-02 

-.23660E-04 

.42403E-04 

-.11091 E-05 

-.78803E-06 

(?) 

0.46 

.14758 E-02 

-.31013E-04 

.5721 6E-04 

-.14187E-05 

-.12807E-05 


0.55 

.17924E-02 

.37951 E-04 

.72349E-04 

-.16251 E-05 

-.16754E-05 

CD 

0.64 

•21287E-02 

-.44881 E-04 

.8S528E-04 

-.1451 8E-05 

-.15477 E-05 

o 

1 

0.73 

.25000E-02 

-.51517E-04 

.1 0952E-03 

-.19970E-05 

-.23488E-05 


0.82 

.29512E-02 

-.57755 E-04 

.1341 IE-03 

-.23460E-0S 

-.28170E-05 


0.91 

.36521 E-02 

-.68062E-04 

.17741 E-03 

-.321 36E-05 

-.40410E-05 


Table 2. Fourier coefficients of the displacement thickness, = .599, 
a = 0°, Re = .48 x 10 , w = 4.789 Hz, Amplitude = 1°. 
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X Location (ft) 

A0 

A1 

B1 

A2 

B2 



.19396E-03 

-.1 581 9E-05 

.53122E-05 

.67836E-07 

.28307E-07 


■ nc 

.36255E-03 

.45777E-05 

-.40789E-05 

-.35501 E-06 

.97161 E-07 

<D 

O 


.54074E-03 

.96275E-05 

-.1 2559E-04 

-.53417E-06 

.91 897E-07 

cd 

: 


.72766E-03 

.14538E-04 

-.21068E-04 

-.66937E-06 

.64907 E-07 

13 

m 

0.46 

.92427 E-03 

.19361 E-04 

-.29917E-04 

-.81875E-06 

•21795E-07 

v/ 

0.55 

.11317E-02 

.24066E-04 

-.39249E-04 

-.931 01 E-06 

-.29324 E-07 

o 

CL 

0.64 

.13532E-02 

.28596E-04 

-.49537E-04 

-.10316E-05 

-.12977E-06 

0 . 

0.73 

.15977E-02 

.32907E-04 

-.61394E-04 

-.12272E-05 

-.22583E-06 

— / 

0.82 

.18904E-02 

.37306 E-04 

-.76396E-04 

-.14735E-05 

-.40425E-06 


0.91 

.23236 E-02 

.43333 E-04 

-.10003 E-03 

-.18755E-05 

-.71813E-06 




.19391 E-03 

.15932E-05 

-.53138E-05 

.54570 E-07 

.87310E-07 

<D 

O 

3 

r 

.36256 E-03 

-.45687E-05 

.40855E-05 

-.32827E-06 

.50857 E-07 

m 

.54082E-03 

-.96154E-05 

.12566E-04 

-.46961 E-06 

-.53037 E-07 

w 

=3 

labile rMI Jlltt 

.72781 E-03 

-.14520E-04 

.21 071 E-04 

-.56848E-06 

-.17724E-06 

CO 

0.46 

.92449 E-03 

-.19332E-04 

.29916E-04 

-.67790E-06 

-.321 97E-06 

<D 

0.55 

.1 1320E-02 

-.24020E-04 

.39238E-04 

-.75995E-06 

-.4781 5E-06 

5 

0.64 

.13536 E-02 

.2851 6 E-04 

.4951 2 E-04 

-.83094E-06 

-.69132E-06 

o 

—i 

0.73 

.15983 E-02 

-.32778 E-04 

.61346E-04 

-.99498E-06 

-.91953E-06 


0.82 

.1891 IE-02 

-.37094 E-04 

.76321 E-04 

-.12185E-05 

-.12506E-05 


0.91 

.23246E-02 

-.43001 E-04 

.99899 E-04 

-.15864E-05 

-.17971 E-05 


Table 3, Fourier coefficients of the momentum thickness, M = *599, 

CO 

a = 0°, Re = .48 x 10^, a> = 4.789 Hz, Amplitude = 1°. 
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Y Location (ft) 

A0 

A1 *1 0E-02 

B1 *1 0E-02 

A2M0E-09 

B2*1 0E-04 

.79905E-05 

.04948 

.06511 

-.02223 

.11204 

.19723 

.19971 E-04 

.12440 

.16135 

-.05247 

.27655 

.49104 

.37933E-04 

.22848 

.27689 

-.07411 

.45449 

.85434 

.64861 E-04 

.33992 

.35805 

-.05323 

.53970 

1.10370 

.10523E-03 

.43214 

.40067 

-.01016 

.53985 

1.20040 

.16575 E-03 

.50287 

.42780 

.03013 

.51757 

1.23020 

.25645 E-03 

.55880 

.44991 

.06339 

.48891 

1.22900 

.39237E-03 

.60619 

.47085 

.09180 

.45228 

1.20170 

.59630E-03 

.64962 

.49201 

.11867 

.40013 

1.14100 

.90094E-03 

.69217 

.51366 

.14719 

.32432 

1.03670 

.13572E-02 

.73587 

.53623 

.17903 

.22760 

.89435 

.20395 E-02 

.78255 

.56236 

.21255 

.14362 

.76938 

.30582 E-02 

.83504 

.59775 

.24339 

.14978 

.78787 

.45759 E-02 

.89799 

.64677 

.27354 

.32806 

1.07040 

.68298E-02 

.97582 

.68939 

.34904 

.64044 

1.50270 

.10161 E-01 

1.05930 

.60300 

.67043 

1.07000 

2.11460 

.1 5052E-01 

1.08810 

.30899 

1.24970 

-.33646 

-.10069 

.22161 E-01 

1.08860 

.29894 

1.26970 

-.50781 

-.33409 

.32340 E-01 

1.08860 

.29850 

1.26970 

-.50740 

-.33452 

.46615E-01 

1.08860 

.29791 

1.26960 

-.50683 

-.33511 

.66057 E-01 

1.08860 

.29714 

1.26940 

-.50608 

-.33586 

.91 504 E-01 

1.08860 

.29622 

1.26930 

-.50656 

-.33494 

.12313E+00 

1.08860 

.29663 

1.26910 

-.53360 

-.29539 


Table 4. Fourier coefficients of the streamwise 

velocity component at oc » 0.64 ft, = .599, 
a = 0°, Re = .48 x 10 , w = 4.789 Hz, Amplitude = 1° . 

7 oo 
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X Location (ft) 

AO 

A1 

B1 

A2 

B2 


0.00667 

.37359 E+00 

.53143E-01 

-.15568E+00 


.10341 E-01 


0.02 

-.16407E+00 

.66513E-01 

-.19742E+00 


.11 953E-01 


0.10 

-.50579E+00 

.34666E-01 

-.12677E+00 


.66401 E-02 


0.20 

-.48913E+00 

.15502E-01 

-.86272E-01 


.36559E-02 

o 

a 

0.30 

-.42419E+00 


-.62910E-01 


.23653E-02 

a 

"C 

0.40 

-.35098E+00 

HBESaiS 

-.47713E-01 


.16619E-02 

3 

CO 

0.50 

* -.27876E+00 

-.46788E-02 

-.37000E-01 


.12288E-02 

ftp* 

0.60 

-.20939 E+00 

-.66793E-02 

-.28956E-01 


.92628E-03 

<D 

CL 

0.70 

-.14071 E+00 

-.75418E-02 

-.22558E-01 


.69629E-03 

CL 

0.80 

-.66016E-01 

-.75031 E-02 

-.17173E-01 


.50170E-03 

-J 

0.90 

.32960E-01 

-.67228 E-02 

-.12298E-01 


.31814E-03 


1.00 

.13490 E+00 

-.65013E-02 

-.11328E-01 


.27956E-03 


2.35 

.25901 E-01 

-.32397 E-03 

-.69127E-05 

.33794E-04 

-.58784E-04 


4.259 

.64130E-02 

-.39253 E-03 

-.14794E-03 

.3609 IE-04 

-.1 1784E-04 





0.00667 

.37234E+00 

-.52208 E-01 

.15568 E+00 

.63 07 6 E-02 

-.90089E-03 


0.02 

-.16566E+00 

-.65501 E-01 

.19744E+00 

.66513E-02 

-.23155E-02 


0.10 

-.50644E+00 

-.34047E-01 

.12677E+00 

.31258E-02 

-.25283 E-02 


0.20 

-.48916E+00 

-.15117E-01 

.86259E-01 

.11 151 E-02 

-.25727 E-02 

O 

o 

0.30 

-.42390E+00 

-.47297E-02 

.62885E-01 

.41 961 E-03 

-.21 672E-02 

CO 

0.40 

-.35051 E+00 

.13454E-02 

.47 677 E-01 

.96358E-04 

-.17663E-02 

=3 

0.50 

-.27820E+00 

.49352E-02 

.36959E-01 

-.56707E-04 

-.14269E-02 

CO 

0.60 

-.20880E+00 

.69319E-02 

.28902E-01 

-.12427E-03 

-.1 1465E-02 

w 

o 

0.70 

-.14015E+00 

.77983E-02 

.22496E-01 

-.14360E-03 

-.91 693E-03 

5 

0.80 

-.65507E-01 

.77729 E-02 

.17100E-01 

-.13505E-03 

-.72658E-03 

o 

—I 

0.90 

.33374E-01 

.70029 E-02 

.12215E-01 

-.99902E-04 

-.56533E-03 


1.00 

.13529E+00 

.6771 1 E+02 

.11 241 E-01 

-.8901 2E-04 

-.53487E-03 


2.35 

.25919E-01 

.38366E-03 

.10067E-03 

.20206E-04 

-.5S207E-04 


4.259 

.64370E-02 

.46165E-03 

.10643E-03 

.1 81 69E-04 

-.221 76E-04 


Table 5. Fourier coefficients of the outer edge 

free stream pressure coefficient, = .599, 
a = 0°, Re = .48 x 10 , w = 4.789 Hz, 
Amplitude = 1 . 
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Figure 1. First Sweep Coordinates 
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Figure 2 . Second Sweep Coordinates 


= 0.599 
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Figure 5 . Three-Dimensional Computational Domain. 




Figure 6. - Three-Dimensional Computational Domain for Nonorthogonal Coordinates 
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Figure 7. Two-Dimensional Section 
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Figure 8. Two-Dimensional Section with Wake 
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Figure 9. 3-D Perspective of Wing Side View 
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Figure 10. 3-D Perspective of Wing Top View 
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Figure 11. 3-D Perspective with Wake 
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Figure 12. Spanwise Coordinates 
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Figure 13 
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